Pharmacophore Based Virtual Screening Approach to Identify Selective PDE4B Inhibitors.

Phosphodiesterase 4 (PDE4) has been established as a promising target in asthma and chronic obstructive pulmonary disease. PDE4B subtype selective inhibitors are known to reduce the dose limiting adverse effect associated with non-selective PDE4B inhibitors. This makes the development of PDE4B subtype selective inhibitors a desirable research goal. To achieve this goal, ligand based pharmacophore modeling approach is employed. Separate pharmacophore hypotheses for PDE4B and PDE4D inhibitors were generated using HypoGen algorithm and 106 PDE4 inhibitors from literature having thiopyrano [3,2-d] Pyrimidines, 2-arylpyrimidines, and triazines skeleton. Suitable training and test sets were created using the molecules as per the guidelines available for HypoGen program. Training set was used for hypothesis development while test set was used for validation purpose. Fisher validation was also used to test the significance of the developed hypothesis. The validated pharmacophore hypotheses for PDE4B and PDE4D inhibitors were used in sequential virtual screening of zinc database of drug like molecules to identify selective PDE4B inhibitors. The hits were screened for their estimated activity and fit value. The top hit was subjected to docking into the active sites of PDE4B and PDE4D to confirm its selectivity for PDE4B. The hits are proposed to be evaluated further using in-vitro assays.


Introduction
Prevalence of inflammatory diseases of respiratory tract i.e. asthma and COPD has increased in recent years, with more than 200 million people affected by it worldwide. Most of the mortality related to these inflammatory disorders occurs in low and low-middle income countries (1).
Phosphodiesterase 4 (PDE4), a predominant family of phosphodiesterase (PDE) enzymes expressed in immune and inflammatory cells, includes three subtypes PDE4A, PDE4B and PDE4D. Inhibition of PDE4 has been shown to suppress a diverse spectrum of inflammatory responses in-vitro and in-vivo (2)(3)(4)(5). More importantly, many PDE4 inhibitors in development are efficacious in animal models of various inflammatory disorders, such as asthma, COPD, psoriasis, inflammatory bowel diseases, and rheumatoid arthritis (3, 6, 7), as well as in clinical trials for asthma and COPD (8-10).
The development of PDE4 inhibitors has been slowed down due to narrow therapeutic window of most of the compounds. A major reason for their poor clinical results is the consequence of dosing limitation caused by side effects such as nausea and emesis (11). Recent findings in PDE4D knockout mice suggest that an inhibitor with PDE4B selectivity should retain many beneficial anti-inflammatory effects without the unwanted side effects (12). The study also established that PDE4D inhibition is responsible for the dose limiting side effects. Some other studies have proven that selective PDE4B inhibitors have potent anti-inflammatory effects in-vitro and in-vivo. Investigation in ferrets also showed significantly less emesis with this compound compared with the non-selective PDE4 inhibitor cilomilast (13). Thus, PDE4B has been established as an extremely attractive target for design of anti-inflammatory drugs, particularly for asthma and COPD.
The highly conserved catalytic domain of PDE4 isozymes makes the design of inhibitors with PDE4 subtype selectivity a challenging task, nevertheless subtype selective PDE4 inhibitors have recently been described (14,15). Only a few amino acids are non-conserved in N-terminal regulatory domain UCR2 (i.e Phe in PDE4D vs Tyr in PDE4B) and C-terminal domain CR3 (i.e Leu in PDE4D vs Gln in PDE4B) (16,17). These minor differences in the regulatory domains have been exploited to design selective PDE4B or PDE4D inhibitors so far (16)(17)(18).
The availability of PDE4B and PDE4D inhibition data for recently reported PDE4 inhibitors allows the development of pharmacophore models of . Pharmacophore models also help in the identification of structural features which can differentiate between the two receptor subtypes. The information obtained can be used for design of more selective and potent PDE4B inhibitors with hitherto new structures. The pharmacophore models of PDE4B and PDE4D inhibitors can be used to screen databases of drug like compounds in a sequential manner to identify novel leads as selective PDE4B inhibitors. Pharmacophore model based virtual screening has proved to be a useful strategy for identification of novel leads in the past (22-32).
In the present study pharmacophore models of both PDE4B and PDE4D inhibitors has been developed and validated. The pharmacophore models were then used for sequential virtual screening to identify novel selective PDE4B inhibitors. The hits were screened for their estimated activity and fit value. Their selectivity for PDE4B was confirmed by docking studies.

Data set
Selective PDE4B inhibitors belonging to thiopyrano[3,2-d] Pyrimidines,2 arylpyrimidines and triazines class reported recently, along-with their PDE4B and PDE4D inhibitory activities, were used for the present study ( Figure 1) (19)(20)(21). The molecular structures and IC 50 of the above series were taken from the original papers. Numbers used in original papers were used to denote molecules belonging to triazine series while numbers used in original papers for molecules belonging to 2-arylpyrimidine and thiopyrano[3,2-d]Pyrimidine series were triazines class reported recently, along-with their PDE4B and PDE4D inhibitory activities, were used for the present study ( Figure 1) (19)(20)(21). The molecular structures and IC50 of the above series were taken from the original papers. Numbers used in original papers were used to denote molecules belonging to triazine series while numbers used in original papers for molecules belonging to 2-arylpyrimidine and thiopyrano[3,2-d]Pyrimidine series were suffixed with a and b respectively.

Pharmacophore generation
Pharmacophore modeling is the most widely used method for identification of essential structural features required for biological activity. In the present study, HypoGen algorithm was applied to 5 Selective PDE4B Inhibitors identification by virtual screening suffixed with a and b respectively.
3D QSAR pharmacophore modeling Pharmacophore generation Pharmacophore modeling is the most widely used method for identification of essential structural features required for biological activity. In the present study, HypoGen algorithm was applied to build the 3D QSAR pharmacophore models for both PDE4B and PDE4D inhibitors using DS V2.0 software (Accelrys Inc., San Diego, CA, USA) (33).
For the study, 75 molecules, with activity values (IC 50 ) between 3.0 nM and 18755 nM were selected as training set, which was used to engender the hypotheses. The training set selected complies with the requirements specified in the literature. To validate the hypothesis, the test set was prepared using the specified requirements. Test set contains 24 molecules having wide range of activity values. Sketch function of DS was used to sketch the two-dimensional (2D) chemical structures of all molecules which were later converted into 3D structures. Maximum of 250 conformations were generated for each molecule using the best conformation model generation method based on CHARMm force field and Poling algorithm (34). Those conformations with energy higher than 20 kcal/mol from the global minimum were rejected. Molecules with their conformational models were then submitted to DS for generating hypotheses.
Automated 3D QSAR pharmacophores were produced by comparing the PDE4B and PDE4D inhibitory activity values of molecules in the training set separately. This helps in identifying the features that are common with the active compounds, but excludes common features for the inactive compounds within conformational allowable regions of space. Selecting the chemical features is one of the most important steps in generating a pharmacophore. While generating hypotheses, HBA (hydrogen bond acceptor), HBD (hydrogen bond donor), and H (hydrophobic), features were selected based on the training set molecules. The number of features allowed in the model were kept in the range 0-5. The 'Uncertainty' values for all the 75 molecules in the training set were set to 2.0, and the default values for other parameters were kept constant. Subsequently, pharmacophore models were computed and the 10 top scoring hypotheses for both PDE4B and PDE4D inhibition were selected for further study. The qualities of the hypotheses were reliant on the fixed cost, null cost, and total cost values (35).
Assessment of pharmacophore quality Quality of the developed pharmacophore was assessed using three different methods. Initially, cross validation was performed by the Fischer's randomization test. Secondly, the prediction of the activity values of the test set was performed. The correlation between the experimental and predicted activities was used to assess predictive ability of the model. All queries were addressed using the ligand pharmacophore mapping protocol.

Virtual Screening
The validated pharmacophore model (Hypo1B and Hypo1D) of PDE4B and PDE4D inhibitors was used as a query in a sequential manner to search the zinc database. Zinc is a comprehensive database of small molecules containing a total of 17,900,742 drug like molecules (36). In the first step ligand pharmacophore mapping module of DS was used along with Hypo1B as the pharmacophore model and zinc database as the database. In the next step, hits mapping to the pharmacophore model Hypo1B were retrieved and hit compounds showing Hypo1B estimated IC 50 less than 20 nM were selected and subsequently subjected to screening using the validated pharmacophore model Hypo1D in the same manner as in the previous step. The hit compounds were chosen that showed Hypo1D estimated fit value less than 4.
Docking studies were used to confirm the selectivity of the hits obtained using the pharmacophore based virtual screening. The most PDE4B selective hit determined by the fit values for Hypo1B and Hypo1D and the most selective ligand from the series used for pharmacophore development i.e. 34b, were docked into the active sites of PDE4B (PDB ID: 4NW7) and PDE4D (PDB ID: 1Y2B). First the protein structures were prepared using the automatic protein preparation module of DS V2.0 software using the default parameters. The structures of the identified hit as well as the standard molecule (34b) were prepared using the prepare ligand module of DS V2.0 software. Docking of the prepared ligands into the active site of the prepared structures of PDE4B and PDE4B was carried out using CDOCKER program available in DS V2.0 software (Accelrys Inc., San Diego, CA, USA) with default parameters (37). The ratio of PDE4B/PDE4D docking scores was used as measure of PDE4B selectivity. The higher is the ratio the greater is the PDE4B selectivity.

Results and Discussion
3D QSAR pharmacophore modeling Pharmacophore generation The top scoring model (Hypo1B) for PDE4B inhibition consist of three HBA which established the highest cost difference (143.378), best correlation coefficient (0.9571), maximum fit value (5.8678) and lowest root mean square (RMS) of 1.86 (Table 1). The results revealed the importance of HBA in PDE4B receptor antagonist activity. The fixed and the null cost values were 236.38 and 509.10, respectively (Table 1). Difference between these two costs (143.378) was greater than 70 bits which showed that the model has over 90% statistical significance. A good pharmacophore model should also have the configuration cost lower than 17, and it was found to be 12.53 for the generated pharmacophore hypotheses. Hypo1B showed correlation coefficient value of 0.9571, demonstrating its good prediction ability.
Top scoring model (Hypo1D) for PDE4D inhibition consists of two HBA and three H with highest cost difference (164.419), best correlation coefficient (0.9563), maximum fit value (8.1515), and lowest root mean square (RMS) of 1.66 (Table 2). As in the case of Hypo1B, HBA was found to be important for PDE4D receptor antagonist activity although there is additional H in this case. Difference between fixed and null costs (164.419) showed that the model has over 90% statistical significance. The configuration cost was also sufficiently low at 12.49. Hypo1D showed correlation coefficient value of 0.9563 (Table 2). Based on statistical parameters Hypo1B and Hypo1D were selected as the best hypothesis for PDE4B and PDE4D inhibition respectively and were employed for further analyses. Figure 2 shows Hypo1B, and Hypo1D chemical features with their geometric parameters while Molecules with highest and lowest activity in the training set aligned to Hypo1B and Hypo1D are shown in Figure 3. The prediction accuracy of both the models was verified using the training set and the activity of each molecule in training set was estimated by regression analysis. The experimental and predicted activities by Hypo1B and Hypo1D for 75 training set molecules are shown in Tables 3 and 4 respectively. Data clearly shows the good agreement between predicted and experimental IC 50 values. Difference between the predicted and experimental values. '+' indicates that the predicted IC 50 is higher than the experimental IC 50 ; '-' indicates that the predicted IC 50 is lower than the experimental IC 50 ; a value of 0 indicates that the predicted IC 50 is equal to the experimental IC 50 .
Close examination of the pharmacophore models Hypo1B and Hypo1D reveals the structural features of an inhibitor which can differentiate well between the two receptors.          The conformation which can allow -COOH at R 3 and hydrophobic groups like halogen atoms in the aromatic ring (Ar) to orient properly for interaction with CR3 will show significant selectivity for PDE4B as compared to PDE4D. This is consistent with the findings described previously in the original papers in which these compounds have been reported (21).

Validation of Hypo1B and Hypo1D
The generated hypotheses were validated using standard methods to check whether the best hypotheses are statistically significant and have considerable predictive ability.

Fischer's randomization method
Fischer's randomization was used to evaluate the statistical significance of the Hypotheses. Validation was done by generating random spreadsheets for training set molecules, which randomly reassigned activity values to every molecule and subsequently generated the hypotheses using the same features and parameters originated for Hypo1B and Hypo1D. All the randomly generated spreadsheets had higher total cost values and lower correlation coefficient values as can be seen clearly from Figure 4. This suggests that Hypo1B and Hypo1D were not generated by chance.

Test set
Test set was prepared using the same protocol as training set and used to determine whether the hypotheses were able to predict the active molecules other than those present in the training set.
The correlation coefficient (r) for the test set given by Hypo1B was 0.8579 (Table 5) while that by Hypo1D was 0.8299 (Table 6). Test set molecules were classified using the same criteria as used for training set molecules. Thus Hypo1B and Hypo1D were able to estimate the PDE4B and PDE4D inhibition activities respectively with reasonable accuracy.

Virtual Screening
Zinc, a comprehensive database of small drug like molecules was used for the sequential virtual screening using the pharmacophore models. Screening of zinc database using the validated pharmacophore model (Hypo1B) of PDE4B inhibitors retrieved a set of 6015 hits, mapping to the pharmacophore model Hypo1B. The hits comprised of some compounds structurally similar to that of the existing PDE4B inhibitors, and some novel scaffolds.
The 397 hit compounds showing Hypo1B estimated IC 50 less than 20 nM were selected and subsequently subjected to screening using the validated pharmacophore model Hypo1D. 5 hit compounds that showed Hypo1 PDE4D estimated fit value less than 4 were identified ( Figure 5). Among the hits ZINC09157416 demonstrated the best PDE4B selectivity based on the hit values (Table 7). ZINC09157416 aligned with Hypo1B and Hypo1D is shown in Figure 6.
The results of docking studies of ZINC09157416 and 34b with PDE4B and PDE4D further confirmed the selectivity of ZINC09157416 for PDE4B over PDE4D (Table 8). randomly generated spreadsheets had higher total cost values and lower correlation coefficient values as can be seen clearly from Figure 4. This suggests that Hypo1B and Hypo1D were not generated by chance.

Test set
Test set was prepared using the same protocol as training set and used to determine whether the hypotheses were able to predict the active molecules other than those present in the training set.
The correlation coefficient (r) for the test set given by Hypo1B was 0.8579 (Table 5) while that by Hypo1D was 0.8299 (Table 6). Test set molecules were classified using the same criteria as used for training set molecules. Thus Hypo1B and Hypo1D were able to estimate the PDE4B and PDE4D inhibition activities respectively with reasonable accuracy. Table 5. Actual and estimated activity of the test set molecules based on the pharmacophore model Hypo1B. Table 5. Actual and estimated activity of the test set molecules based on the pharmacophore model Hypo1B. showing Hypo1B estimated IC 50 less than 20 nM were selected and subsequently subjected to screening using the validated pharmacophore model Hypo1D. 5 hit compounds that showed Hypo1 PDE4D estimated fit value less than 4 were identified ( Figure 5). Among the hits ZINC09157416 demonstrated the best PDE4B selectivity based on the hit values (Table 7).

21
ZINC09157416 aligned with Hypo1B and Hypo1D is shown in Figure 6. The results of docking studies of ZINC09157416 and 34b with PDE4B and PDE4D further confirmed the selectivity of ZINC09157416 for PDE4B over PDE4D (Table 8).

Conclusions
Ligand-based pharmacophore models for a diverse class of PDE4B and PDE4D inhibitors were developed. The best pharmacophore models Hypo1B and Hypo1D were validated using different methods to evaluate their predictive power over the diverse test set compounds. Hydrogen bond acceptors were identified to be mainly responsible for PDE4B inhibition while both hydrogen bond acceptors as well as hydrophobic groups were found to be responsible for PDE4D inhibition. The highly predictive pharmacophore hypotheses were further used in sequential virtual screening for identification of selective PDE4B inhibitors.
Zinc drug like database was used in virtual screening. The hits from the virtual screening were filtered based on the estimated activity and fit value. Five molecules with different backbones were identified as final hits. The most selective hit molecule ZINC09157416 exhibited better selectivity for PDE4B than the standard compound 34b in the docking studies. The activity of the hit compound has not been reported in the literature as we explored by PubChem and SciFinder Scholar search tools. Thus, the sequential virtual screening strategy using 3D QSAR pharmacophores for PDE4B and PDE4D inhibitors proved to be an effective strategy to identify novel selective PDE4B inhibitors.